Photoacoustic image reconstruction from 
few-detector and limited-angle data 

Lei Yao and Huabei Jiang* 

Department of Biomedical Engineering, University of Florida, Gainesville, FL 32611, USA 

*hjiang@hme. ufl. edit 

Abstract: Photoacoustic tomography (PAT) is an emerging non-invasive 
imaging technique with great potential for a wide range of biomedical 
imaging applications. However, the conventional PAT reconstruction 
algorithms often provide distorted images with strong artifacts in cases 
when the signals are collected from few measurements or over an aperture 
that does not enclose the object. In this work, we present a total-variation- 
minimization (TVM) enhanced iterative reconstruction algorithm that can 
provide excellent photoacoustic image reconstruction from few-detector and 
limited-angle data. The enhancement is confirmed and evaluated using 
several phantom experiments. 
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1. Introduction 

Photoacoustic tomography (PAT) is an emerging non-invasive imaging technique that 
combines the merits of high optical contrast and high ultrasound resolution in a single 
modality [1-3]. In PAT, the Helmholtz-like photoacoustic wave equation has been commonly 
used as an accurate model for describing laser-induced acoustic wave propagation in tissue. 
While a variety of analytic PAT reconstruction algorithms have been developed [4—6], the 
finite element method (FEM) based approach appears to be particularly powerful in this 
regard, because it can provide quantitative imaging capability by recovering optical absorption 
coefficient [7,8], eliminate the assumption of homogeneous acoustic medium needed in 
analytical methods, accommodate object boundary irregularity and allow appropriate 
boundary conditions implementations. 

A practical need exists for reconstruction of photoacoustic images from few 
measurements, as this can greatly reduce the required scanning time and the number of 
ultrasound sensors placed near or on the boundary of an object to receive the laser-induced 
acoustic signals. In addition, in many practical implementations of PAT the photoacoustic 
signals are recorded over an aperture that does not enclose the object, which results in a 
limited-angle tomographic reconstruction problem. In such cases, the existing reconstruction 
algorithms, which are based on the least-squares criteria (i.e., the regularized Newton method) 
[9,10], often generate distorted images with severe artifacts. Total-variation-minimization 
(TVM), on the other hand, has proved to be a powerful tool for limited-data image 
reconstruction in diffraction tomography and computed tomography (CT) [11,12]. TVM based 
PAT algorithms have also been implemented and tested using numerically simulated limited- 
view data [13-16]. Here we describe an unconstrained TVM FEM-based iterative 
reconstruction algorithm, and for the first time present experimental evidence that the TVM- 
based algorithm offers excellent photoacoustic image reconstruction from few-detector and 
limited-angle data. 

This paper is organized as follows. In Section 2, the FEM based PAT reconstruction 
algorithm and the unconstrained TVM scheme are reviewed briefly. The experimental 
validation of our TVM FEM-based algorithm is presented in Section 3. The conclusions are 
made in Section 4. 

2. Method 

We first briefly describe the existing FEM-based photoacoustic reconstruction algorithm 
detailed elsewhere [17]. The time-domain photoacoustic wave equation in tissue can be 
described as follows: 

P[ ' ' vl dt 2 C p dt ' ( ) 

where p is the pressure wave; v 0 is the speed of acoustic wave in the medium; /? is the thermal 
expansion coefficient; C p is the specific heat; <D is the absorbed energy density; 
J (/) = S(t - 1 0 ) is assumed in our study. 
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To form an image from a presumably uniform initial guess of the absorbed energy density 
distribution, we need a method of updating <D, the desired quantity from its starting value. This 
update is accomplished through the least-squares minimization of the following functional: 

M , 

np.°)=i l (p o j-p e j). <2> 

where p° and p c j are the measured and computed acoustic field data for i = \,2,---M 

boundary locations. Using the regularized Newton method, we obtained the following matrix 
equation for updating <D: 

(3 T 3 + M)A Z = 3 T (p°-p c ), (3) 

where p° =(p° , p",-" Pm) an d p '=(K , pW" Pm) > 4£ is tne update vector for the 
absorbed optical energy density; 3 is the Jacobian matrix formed by dp/d<& at the boundary 

measurement sites; I is the regularization parameter determined by combined Marquardt and 
Tikhonov regularization schemes; and I is the identity matrix. 

We now incorporate the total variation of O as a penalty term by defining a new functional 
[18]: 

F(p,0) = F(p,0) + L(0). (4) 

Here i( < 5) = j-\ja>l ft^f + S 2 dxdy is the penalty term, and co^ and S are typically positive 

parameters that need to be determined numerically. The minimization of Eq. (4) can be 

realized by the differentiation of F with respect to each nodal parameter that constitutes the O 
distribution and by setting each of the resulting expression to zero, leading to the following 
system of equations: 




where V t = 3L/3<D ; . Again, through the regularized Newton method, the following matrix 
equation for TVM-based inversion is obtained [18]: 

(3 T 3 + R + M)Az = 3 T (p° -p c )-V, (6) 

where V is formed by 3L/3<D and R is formed by 3V/3<D . At this point, our solution 

procedure for solving Eq. (6) and the regularization parameter selection are identical to those 
described previously. Hence, it becomes clear that the only additions to our new algorithm 
result from the assembly of matrix R and the construction of column vector V. 

3. Results 

In this section our TVM enhanced reconstruction algorithm is tested and evaluated using 
phantom experimental data. For comparative purposes, reconstruction results without the 
TVM enhancement are also presented. 

The experimental setup used for collecting the phantom data was a pulsed ND: YAG laser 
based single transducer (1MHz) scanning system, which was described in detail elsewhere [8]. 
Three phantom experiments were conducted. In the first two experiments, we embedded one 
or two objects with a size ranging from 3 to 0.5 mm in a 50 mm-diameter solid cylindrical 
phantom. The absorption coefficient of the background phantom was 0.01 mm -1 , while the 
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absorption coefficient of the target(s) was 0.03 mm -1 . In the last experiment, we used a single- 
target-containing phantom, aiming to test the capability of detecting target having low optical 
contrasts relative to the background phantom. In this case, the target had an absorption 
coefficient of 0.015 mm -1 . The reduced scattering coefficients of the background phantom and 
targets were 1.0 and 3.0 mm -1 for the first two experiments, and 1.0 and 2.0 mm -1 for the last 
experiment. In the reconstructions, we used a dual-meshing scheme [19] where the fine mesh 
used for the forward calculation consisted of 5977 nodes and 11712 elements, while the 
coarse mesh used for the inverse calculation had 1525 nodes and 2928 elements. All the 
images obtained from the method without the TVM are the results of three iterations, while 
those obtained from the TVM-enhanced method are the results of twenty or more iterations. 
Parallel code was used to speed up these calculations on Beowulf clusters with 8 CPUs. The 
computational speed can be further increased if clusters with more than 8 CPUs are used. As 
mentioned above, the parameters and S were determined through numerical 

experimentation. From our experience, a constant value of S = 0.001 was sufficient for the 
current experimental studies, while the value of a> m is related to the signal-to-noise ratio of the 

measurements. For our experimental cases, co^ =1.0 and 8 = 0.001 for the first case, and 
= 0.5 and S = 0.001 for the second and third cases were used. 

The reconstruction results based on few-detector data from the three experimental cases 
are shown in Figs. 1, 2 and 3, respectively. In each figure, the images in the top and bottom 
rows, respectively, present the recovered absorbed energy density images using the algorithm 
without and with the TVM where 120, 60, 30, and 15 detectors were equally distributed along 
the surface of the circular background region. As expected, the conventional algorithm 
without the TVM provides high quality images only when the number of detectors was 
relatively sufficient (Figs. la,b, Figs. 2a,b, and Figs. 3a,b). The images reconstructed from 
few-detector data by the conventional algorithm without the TVM contained severe artifacts 
and distortions (Figs. lc,d, Figs. 2c,d, and Figs. 3c,d). Considerably enhanced images are 
achieved using the method with the TVM, especially when the number of the detectors is 
reduced to 30 or 15 (Figs. lg,h, Figs. 2g,h, and Figs. 3g,h). We also note that the 
computational efficiency of our TVM based algorithm for recovering a small absorber (Fig. 2) 
and a larger absorber (Fig. 3) is similar when the number of the detectors is insufficient. 

Images reconstructed from the limited-angle data for the first case using the two methods 
are displayed in Fig. 4 where the top and bottom rows, respectively, show the recovered 
absorbed energy density images using the algorithm without and with the TVM when 120 
detectors over 360°, 60 detectors over 180°, and 30 detectors over 90°, were equally 
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Fig. 1. Reconstructed photoacoustic images based on few-detector data for case 1. (a), 120 
detectors, without TVM. (b), 60 detectors, without TVM. (c), 30 detectors, without TVM. (d), 
15 detectors, without TVM. (e), 120 detectors, with TVM. (f), 60 detectors, with TVM. (g), 30 
detectors, with TVM. (h), 15 detectors, with TVM. 
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Fig. 2. Reconstructed photoacoustic images based on few-detector data for case 2. (a), 120 
detectors, without TVM. (b), 60 detectors, without TVM. (c), 30 detectors, without TVM. (d), 
15 detectors, without TVM. (e), 120 detectors, with TVM. (f), 60 detectors, with TVM. (g), 30 
detectors, with TVM. (h), 15 detectors, with TVM. 
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Fig. 3. Reconstructed photoacoustic images based on few-detector data for case 3. (a), 120 
detectors, without TVM. (b), 60 detectors, without TVM. (c), 30 detectors, without TVM. (d), 
15 detectors, without TVM. (e), 120 detectors, with TVM. (f), 60 detectors, with TVM. (g), 30 
detectors, with TVM. (h), 15 detectors, with TVM. 
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Fig. 4. Reconstructed photoacoustic images based on limited-angle data for case 1. (a), 120 
detectors over 360°, without TVM. (b), 60 detectors over 180°, without TVM. (c), 30 detectors 
over 90°, without TVM. (d), 120 detectors over 360°, with TVM. (e), 60 detectors over 180°, 
with TVM. (f), 30 detectors over 90°, with TVM. 

distributed along the surface of the circular background region. Again, strong artifacts and 
distortions exist in the images recovered with the method without the TVM when the data 
collected is angle-limited (Figs. 4b,c), while considerably improved images are reconstructed 
using the TVM enhanced algorithm (Figs. 4e,f). 
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In order to evaluate quantitatively the reconstruction quality using the method with and 
without TVM enhancement, we use the following universal quality index (UQI) [20]: 



UQl\i\f°\ 



2Cov f\f° 



2/'/° 



(7) 



Here the image f can be interpreted as vectors of size N: f = (/j, f 2 ,-"f N ) T , where N denotes 
the number of image data acquired from the FEM based algorithm, f' is the image means, 
a J is the variances, and Covif ,f 0 } is the covariances over the whole image domain, where j 



= 0 and 1 . UQI measures the image similarity between the reconstructed ( f 1 ) and reference 
( f 0 ) images, and its value ranges between 0 and 1 . The value of the UQI is closer to 1 when 
the reconstructed image is more similar to the reference image. 

We calculated the UQIs for the three cases presented above, and the results are given in 
Fig. 5 where Figs. 5a and 5b show the results from the few-detector data based on the method 
without and with the TVM, respectively, while Fig. 5c presents the results from the limited- 
angle data for the first case based on the two methods. Here the image data recovered by 120 
detectors was regarded as the reference one. The computed UQIs shown in Fig. 5 confirm the 
observation that the TVM enhanced PAT algorithm provides significantly better image quality 
compared to the conventional method without the TVM. 





Fig. 5. UQIs calculated from the recovered images for the three cases based on few-detector 
data without TVM (a) and with TVM (b), and for case 1 based on limited angle data with and 
without TVM. 

4. Conclusions 

In this work, we have implemented and evaluated an unconstrained, total-variation- 
minimization method for time-domain FEM based PAT. This study has confirmed that in the 
situation that the data is collected from few measurements or over an aperture that does not 
enclose the object, the developed TVM based PAT algorithm provides considerably improved 
image reconstruction compared to the conventional methods. The application of this TVM 
enhanced reconstruction algorithm will significantly reduce the number of ultrasound 
transducers and scanning time needed for high quality photoacoustic image reconstruction. 
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